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We study the spherical, top-hat collapse model for a mixed dark matter model including cold 
dark matter (CDM) and massive neutrinos of mass scales ranging from m„ ~ 0.05 to a few 0.1 eV, 
the range of lower- and upper-bounds implied from the neutrino oscillation experiments and the 
cosmological constraints. To develop this model, we properly take into account relative differences 
1 between the density perturbation amplitudes of different components (radiation, baryon, CDM and 

neutrinos) around the top-hat CDM overdensity region assuming the adiabatic initial conditions. 
■ Furthermore, we solve the linearized Boltzmann hierarchy equations to obtain time evolution of 

" the lineariezed neutrino perturbations, yet including the effect of nonlinear gravitational potential 

due to the nonlinear CDM and baryon overdensities in the late stage. We find that the presence 
of massive neutrinos slows down the collapse of CDM (plus baryon) overdensity, however, that the 
neutrinos cannot fully catch up with the the nonlinear CDM perturbation due to its large free- 
streaming velocity for the ranges of neutrino masses and halo masses we consider. We find that, 
just like CDM models, the collapse time of CDM overdensity is well monitored by the linear-theory 
extrapolated overdensity of CDM plus baryon perturbation, smoothed with a given halo mass scale, 
if taking into account the suppression effect of the massive neutrinos on the linear growth rate. 
Using these findings, we argue that the presence of massive neutrinos of mass scales 0.05 or 0.1 eV 
may cause a significant decrease in the abundance of massive halos compared to the model without 
' the massive neutrinos; e.g., by 25% or factor 2, respectively, for halos with 1O 15 M0 and at z = 1. 
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Galaxy clusters are the most massive, gravitationally bound objects in the universe, therefore, their abundance and 
its redshift-evolution are very sensitive to cosmology including the nature of dark energy [l|, Q as well as the primordial 
non-Gaussianity [H, Hj]. There are various ongoing and planned surveys capable of finding many clusters under a 
homogeneous and well-controlled/calibrated selection; e.g., optical surveys such as the Sloan Digital Sky Survey 
(SDSS) Q, Subaru Hyper SuprimeCam Survey [6] 64] and the arcminute-resolution cosmic microwave background 
(CMB) experiments, Atacama Cosmology Telescopep} and South Pole Telescope @, with which clusters can be 
found via the Sunyaev-Zel'dovich (SZ) effect. These cluster catalogs can be used to derive stringent constraints on 
cosmological parameters or more generally test the paradigm of cold dark matter dominated structure formation 
^h! scenario @,i,[Ioj]. 

t-H ■ The neutrino oscillation experiments have revealed that the standard three-flavor neutrinos have non-zero masses, 
implying that the Big-Bang relic neutrinos contribute to the present-day mean matter density by at least about 0.4 

• i-H ■ and 0.8% if the neutrinos follow the normal mass and inverted mass hierarchies corresponding to the lower bounds 
on the sum of neutrino masses, 0.05 and 0.1 eV, respectively [e.g.,|Tl|. While the absolute neutrino mass scale is still 
unknown, large-scale structure probes provide a powerful means of constraining the neutrino mass [t| [H, EH • In fact 
the cosmological probes have put currently the most stringent upper bound on the sum of neutrino masses; m^tot < 0.2 
- 0.8 eV (95% C.L.), the different bounds for different probes that employ the different level of assumptions on 
nonlinear structure formation [l|, [Til . Il5| . The ongoing and upcoming cosmological surveys promise to further tighten 
the neutrino mass constraint and have the potential to detect the absolute mass scale, rather than the upper bound, 



if systematic errors are well under control. See [16| for the current status of the cosmological neutrino constraints and 
the future prospect. 

However, the previous cluster cosmology experiments rest on the use of halo mass function calibrated based on 
N-body simulations that ignore the effect of massive neutrinos (e.g., [13, lUjl)- Although there are several attempts to 
simulate nonlinear structure formation in a mixed dark matter (CDM plus massive neutrinos) model (see [l9l |20] | for 
the pioneer work and [2lTl24| for the recently revisited attempts), it is still very difficult to accurately simulate the 
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structure formation especially for the neutrino mass scales of a few 0.1 eV or lighter, because such light mass-scale 
neutrinos have a fast free-streaming motion, larger than the gravity-induced bulk motion, at relevant redshift and it 
is difficult to represent the (perturbed) Fermi-Dirac distribution with a finite number of N-body particles at every 
spatial position [see I22I . for an attempt on the grid- based simulation of neutrinos] . There are also several attempts 
[25l428| aimed at developing the perturbation theory based approach to analytically model the nonlinear structure 
formation for a MDM model by extending the linear perturbation theory [29j. However, the perturbation theory 
based model is only valid up to the quasi nonlinear regime and breaks down for the nonlinear regime relevant for halo 
formation. 

Therefore, the purpose of this paper is to develop a top-hat spherical collapse model for a MDM model fully 
taking into account the effects of multi-component system (radiation, baryon CDM and neutrinos) [jjjj HH for the 
pioneer work on the spherical collapse model for a CDM model] [also see [32], [33| . There are several key ingredients to 
include in developing this model. Firstly, we carefully account for differences in the density perturbation amplitudes 
of each component around the top-hat CDM overdensity region assuming the adiabatic initial conditions as predicted 
from standard inflation scenario. Secondly, to achieve the desired accuracy, we solve time evolution of the density 
perturbations from the deeply radiation-dominated regime to the present time or from the sufficiently linear regime 
to the non-linear regime, where the initial density perturbations are on super- horizon scales [34] . Hence, we properly 
take into account the transition of perturbations from the super- to sub-horizon scales. Thirdly, to study the neutrino 
perturbations around the top-hat CDM overdensity, which cannot be treated as a fluid, we properly solve the linearized 
Boltzmann hierarchy equations [291 ] taking into account the nonlinear gravitational potential well due to the nonlinear 
CDM and baryon perturbations. 

The spherical collapse model is an approximated method for studying the nonlinear dynamics due to the unrealistic 
symmetry assumed. Nevertheless, this gives a very useful tool for studying various effects on nonlinear structure 
formation; the spatial curvature or the cosmological constant [35l. |36| . time- varying and/or clustered dark energy [37i — 
|42| . the modified dark matter scenario (43| , the modified gravity scenario [44j - l46j , and the effect of baryon perturbation 
34]. As for the effect of massive neutrinos, we know that the effect is small given the current upper bounds on the 
sum of neutrino masses, ;$ a few 0.1 eV (at most a few % contribution to the matter density). Hence we expect that, 
once the nonlinear dynamics is realized for a MDM model, we can perturbatively include the effect of massive neutrino 
on the CDM simulation based predictions by slightly modifying the model ingredients. For example, the halo mass 
function is given as a function of the peak height v = S c /a(M, z), where cr(M, z) is the rms linear mass fluctuations 
of halo mass scale M at redshift z and S c is the linear-theory extrapolated critical density as indeed motivated by the 
top-hat spherical collapse model [43] ■ Given these facts, we may be able to infer the effect of massive neutrinos on 
the halo mass function once the effects on 5 C and a(M, z) are realized. Thus, along this approach, we will also discuss 
the impact of massive neutrinos on the abundance of massive halos. 

Throughout this paper, we employ, as our fiducial cosmological model, a flat A-dominated CDM model that is 
consistent with the WMAP 7-year result (48[; the present-day density parameters of matter and baryon are flmoh 2 = 
0.1334 and ftboh 2 = 0.0226, respectively; the dimension-less Hubble parameter is h — 0.71; the spectral tile and 
normalization parameter of the primordial power spectrum are n s = 0.963, and A s = 2.43 x 10~ 9 , respectively. In 
most parts of our paper, when adding massive neutrinos to the fiducial cosmological model, we vary the CDM density 
parameter fi c o by fixing the total dark matter density to Vl c oh 2 + £l V Qh 2 = 0.1108, where the energy density parameter 
of massive neutrino is specified by the sum of neutrino masses as il^oh 2 = m^tot/94.1 eV [ll|. We assume the three 
neutrino species, and assume one species among them is massive for simplicity. 

II. METHODOLOGY 

In this section, we develop a method for solving a spherical top-hat collapse of CDM overdensity in a multi- 
component system, which consists of radiation (R), baryon (b), cold dark matter (c) and massive neutrinos {v). 

A. Evolutionary equation of top-hat overdensity on superhorizon scale 

To achieve a sufficient accuracy of the spherical top-hat collapse, we set up the initial conditions of perturbations in 
a deeply radiation dominated era, Zj = 10 7 in our case. The primary reason of this high initial redshift is to reproduce 
the results in the earlier work [3J] at the limit of massless neutrino case (m u , t ot — > 0), where 34] studied the effect 
of baryon perturbation (without massive neutrinos) on halos relevant for first stars at redshift z ^ 30. Note that, for 
cluster-scale halo formation, we can start from the later initial redshift such as z ~ 10 4 -10 5 , but here we use the initial 
redshift Zi ~ 10 7 in order to retain a broader coverage of our model validity. In such a radiation dominated era, the 
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perturbations of interest are on superhorizon scales. To avoid gauge ambiguities that may arise in dealing with such 
superhorizon-scale perturbations, we employ the following approach. 

Since we are interested in a spherical top-hat overdensity region, time evolution of such a top-hat region can be 
most readily described by a perturbed Friedman universe (e.g., 33]). Here we meant by "perturbed" that the top- hat 
region is described by a Friedman universe with a small positive curvature that corresponds to the top-hat density 
perturbation, where the top-hat region is embedded in the background Friedman universe that has a flat geometry. 
Thus the equation of motion for the boundary radius of top-hat overdensity region is given as 
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[pcb(l + Scb) + pr(1 + S R )} - -jg 



(1) 



where R is the radius of the top-hat overdensity region of matter (CDM and baryon) and radiation, p c b and pR 
are the mean energy densities of matter (CDM plus baryon) and radiation, and 5 c b and Sr are their overdensities 
(e.g., defined as 5r = Pr/Pr — 1), respectively. The dot notation ' denotes the time derivative, and the constant k 
is the effective curvature parameter (fc > 0) which is given in terms of the initial overdensity (see below). In this 
regime, massive neutrinos with mass scales we are interested in are still relativistic and contribute the radiation. Note 
S c b = {p~c/pm)S c + (p~b/pAi)8b — 5 C because of 5 C = 5b on superhorizon scales for adiabatic initial conditions, where pu 
is the mean density of total matter (pM = pc + Pb or pM = pc + Pb + Pv when the massive neutrino is relativistic or 
non-relativistic, respectively). Also note that dark energy contribution is negligible in a radiation dominated regime. 
The scale factor for the background universe obeys 
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We use the linear theory predictions to determine the initial conditions of perturbations. The spherical top-hat 
collapse model is equivalent to the case that the perturbations are solved under the synchronous gauge condition. In 
this case, the superhorizon-scale perturbations grow as 5 c b ^ o? [29L l34l]. Assuming this growing mode and using the 
mass conservation R 3 p c b(l + 5 c b) = constant yield the initial condition for the velocity R/R: 



2 Orb i 

3 1 + 5 c b, 



(3) 



where Hi = H(ti) and 5 c b,i = 5 c b(ti). By inserting Eq. Q into Eq. (TT]) at the initial time ti we can express the 
effective curvature parameter k in terms of the initial overdensity 5 cb ^ as 
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where we have used the adiabatic initial condition to re-express radiation perturbation in terms of matter perturbation 
as 5r ~ (3/4)£ c fc. Hence Eq. (fTJ can be re- written as 
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The relation between R and <5 follows from the mass conservation; (1 + 5 c b) = (R/Ri) 3 (1 + 5 c b t i)(a/ai) 3 . Hence, 
given the initial conditions on Ri (or 5 c b t i) and Eq. ([5]) can be solved numerically to obtain time evolution of the 
top-hat overdensity 5 c b until the top-hat region enters into the horizon. 

Eq. is an exact equation that can be applied even if the perturbation amplitude is large (unrealistic though) . In 
the radiation dominated regime, the linear theory gives a good approximation. By linearizing Eq. ([5]) we can derive 
the differential equation which governs time evolution of the linear perturbation <5 L : 



H5^ b = -4ttG 
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where we have used the fact pr^ S> pu,i at the initial redshift. Again note 5 c b 
regime. 



5^ b to a good approximation in this 
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B. Evolutionary equations of perturbations on sub- horizon scales 

When the overdensity region enters into the horizon, perturbations of different components evolve in different ways; 
the CDM overdensity continues to grow, and baryon and neutrinos cannot grow together with CDM. We below 
describe our treatments of each component's dynamics on subhorizon scales after the horizon crossing. 

1 . CDM perturbation 

CDM plays a major role in the spherical collapse model. When the top-hat overdensity region enters into the 
horizon, we use the following equation, obtained in the Newtonian gauge, in order to solve the dynamics up to the 
nonlinear collapse of spherical CDM overdensity region: 

R c (t) 47rG r _ - , GSM(< R c ;t) 

where R c (t) is the radius of the top-hat region of CDM overdensity, and ptot and Ptot denote the mean energy and 
pressure densities, which determine the cosmic expansion history over the range of radiation, matter and dark energy 
dominated eras. Note i? c (£ ntcr) = -Reenter) at the horizon crossing, where R(t enter ) is the radius of the initial top-hat 
overdensity region discussed in the preceding section. The quantity SM is the mass fluctuation within the CDM- 
overdensity sphere, and includes contributions from CDM, baryon and neutrino perturbations: SM = SM c + SMt, + SM v . 
Note that we ignore perturbations of dark energy in this paper [see H^, IU for the spherical collapse model with dark 
energy perturbations, but without massive neutrinos]. 

For CDM perturbation, the mass conservation within the top-hat region holds: 

AttR 3 c _ ^R%i .enter _ , x , ^ R l.i. , Q , 
M c = g Pc\t)[l + c (t)\ = — Pe^cntcrU + Oc,*, enter J — ^~ Pc.i, (») 

where S c (t) is the overdensity at time t, and the quantities with subscript " c ,i,enter" denote their quantities at the 
horizon crossing. For the perturbations of interest, the horizon crossing is earlier than the decoupling epoch, where 
the perturbations are well in the linear regime. In the equation above we have used S cA <C 1 at the initial redshift, 

so used the fact p c ,i(l + S c ,i) — Pc,i- We determine 8 c ,i, enter, Rc,i,entei and R c ,i, enter by matching the values to those 
from Eq. ([5]) at the horizon crossing. We again stress that, by matching these initial conditions from superhorizon 
to subhorizon scales, we can achieve a sufficiently accurate setup of the initial conditions needed for the nonlinear 
spherical collapse dynamics, even for a high collapse-redshift such as z ^ 30. 

To solve Eq. |[7J, we need to specify the mass fluctuations of baryon and massive neutrinos, SM^t) and 6M u (t), 
within the spherical region of radius R c (t) at each time step. However, unlike CDM, the mass fluctuations are not 
conserved within the CDM top-hat region, because baryon is dragged out of the CDM potential well due to the tight 
coupling with photons before the decoupling epoch, and neutrinos are free-streaming out of the CDM potential well 
due to their large thermal velocities [65[ . In fact the linear perturbation theory gives S c 3> o~& S v at the decoupling 
epoch. Hence, for simplicity, we assume 6b = 5 V = during epochs after the horizon crossing to the decoupling epoch. 
One can then find that Eq. (0) gives growing modes of S c oc In a and oc a in the radiation and matter dominated eras 
in the linear regime (i.e. 5 C <C 1), as expected from the linear perturbation theory. 

After the decoupling epoch (z ~ 10 3 ), baryon becomes cold, and follows the mass conservation. We will below 
describe our treatments of baryon and neutrino perturbations in subsequent subsections. 

Linearizing Eq. ([7]) yields the following equation to describe the evolution of linear CDM perturbation: 

5 C L + 2HS^ - 4nG [p c S^ + p b 5l{< i? c L ) + #c)] - 0, (9) 

where 5^ is the linear CDM overdensity, and 8\ (< R c ) and 8^(< R c ) are the linear density perturbations averaged 
within the sphere of radius R^. Here the comoving radius of R% is set to the same in the linear theory: R^/a = R\ A j aj. 

The initial conditions for 8\ i and S^ A are set by matching to Eq. ([6]) at the horizon crossing. Before the decoupling- 
epoch and until sufficiently higher redshift before halo formation, the equation above gives a good approximation to 
Eq. 0. 

2. Baryon perturbation 



Next let us consider evolution of baryon perturbation in a CDM top-hat overdensity region. 
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After the horizon-crossing of the spherical top-hat overdensity region until the decoupling epoch (zdcc), baryon is 
tightly coupled to photon and the baryon density perturbation cannot grow. More exactly, the baryon-photon fluid 
oscillates according to the acoustic sound wave in the CDM potential well - the baryon acoustic oscillations (BAO) 
[49j . The characteristic scale of this clustering is about 150 Mpc for our fiducial cosmological model. Therefore, even 
if the baryon perturbation initially had a top-hat overdensity profile as in CDM, the baryon perturbation becomes 
increasingly spatially-extended than the CDM top-hat region as time goes by until Zdcc- The linear perturbation 
theory predicts that the baryon density perturbation becomes much smaller in the amplitude than the CDM density 
perturbation at the decoupling epoch for length scales of interest; 5b <C 5 C at z^ ec for scales of interest. 

Therefore we simply assume that the baryon density perturbation averaged within the CDM top-hat region is 
5},(< R c ) =0 during epochs from the horizon-crossing to the decoupling epoch; z en t cr > z > z^ cc , where z^ cc is specified 
once a background cosmological model is given, e.g., from the CAMB code [50]. We have checked that the spherical 
collapse of CDM overdensity is not changed even if we instead use the different assumption; Sb(< R c ) — h{< Rc', Center), 
where 5b(< R c ; z cn tcr) is the baryon overdensity when the CDM top-hat region had the horizon-crossing. 

After the decoupling epoch baryon becomes a cold component, and can cluster together with CDM perturbation. 
We can thus use the spherical collapse equation to solve nonlinear evolution of baryon perturbation in the CDM 
top-hat region just like the CDM case (Eq. [7]): 

R b (t) 4ttG f _ - , G5M(<R b -t) 

Here Rb(t) is the radius of baryon overdensity region, which is chosen from the radius satisfying Rb = R c at Zd cc - After 
the decoupling epoch, the mass conservation holds: Rb(t) 3 pb(t) [1 + 5b(< Rb)] — constant. The initial conditions are 
set to R b — HRb at Zd cc , which corresponds to 6b — 0, motivated by the fact that baryon-photon coupling prevents 
baryon perturbation from growing before the decoupling epoch. Thus the baryon velocity perturbation is different 
from that of CDM perturbation, given as Rb > R c , so the time evolution of baryon radius Rb(t) differs from the CDM 
radius R c (t). As times goes by, the baryon perturbation eventually catches up with the CDM perturbation as we will 
explicitly show below. 

3. Neutrino perturbation 

The neutrino perturbation on sub-horizon scales cannot be captured by the CDM potential well due to the large 
free-streaming velocity. As a result, the neutrino perturbation around the CDM overdensity extends out to radius 
comparable with the free-streaming scale that is given as Af s ~ a~ 1 H~ 1 (j VM (z) in the units of comoving scale; here 
a v>v is the velocity dispersion of the Fermi-Dirac distributed neutrinos (see Appendix A in [111 ] for the definition). 
For example, for neutrinos of mass scale m„ ~ 0.1 eV, A ~ 20 h^ 1 Mpc at z = 0. To model these physical processes 
we use the modified CAMB code to solve time evolution of neutrino clustering, where we properly take into account 
the nonlinear gravitational potential due to nonlinear CDM and baryon perturbations in the late stage (see [25| for 
the similar approach to solving the evolution of mildly nonlinear perturbations) . 

To be more precise, we solve the linearized Boltzmann equation that governs time evolution of the neutrino distri- 
bution function f(x\q, hj, r) = /o(e)[l + ^(x l , q, hj,r)]. Here r is the conformal time, q and fij denote the comoving 
momentum and its direction, e = (q 2 + a 2 rn 2 ,) is the proper energy times scale factor a(t), fo is the background 
distribution function (the Fermi-Dirac distribution) and \1/ is the perturbed distribution function. The Boltzmann 
equation in an expanding universe can be reduced to the following hierarchical equations in Fourier space |29j : 

*° = h , 

e 6 a m q 

qk , ek d In f n , „ 

*1 = ^- *o - 2* 2 - -V-r^ , (11) 
6e oq ding 

*i = (2 / + fc 1)£ [Mi-i -V+ (I > 2) , 

where ip and <f> are the metric perturbations in the Newtonian gauge, and the perturbed distribution function * is 
expanded in terms of the Legendre polynomials as 



h, q, t) = ^{-i) e (2£ + l)*,(fc, q, r)P e (k ■ n). (12) 
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Our main interest is to study the impact of massive neutrinos on the spherical collapse of CDM overdensity in the 
matter or dark energy dominated era. At redshifts after the decoupling epoch, the difference between metric pertur- 
bations ifi and (f), which arises from anisotropic stress, is negligible: ip = <f>. In this case, the potential perturbation is 
given by the Poisson equation as 

- k 2 ip{k, t) = AnGa 2 ^ p t (T)6 l (k, r), (13) 

i— b,c,i/ 

where Si(k) denotes the Fourier-transformed coefficients of the i-th component (i —c, b, v). We insert the nonlinear 
top-hat overdensities of CDM and baryon into the Poisson equation above to compute the potential including the 
nonlinear contribution. Note that we take the center of the top-hat region as the coordinate center, which makes the 
potential dependent only on the length of wavevector fc; ip(k) = ip(k). The corresponding potential for the linear- 
theory extrapolated density perturbations is computed from — k 2 ip h (k,T) — AnGa 2 J2i=b c v Pi( r )^H^> t) . Since 
we have assumed that CDM and baryon overdensities take spherical top-hat profiles, whose radii are R c and Rb, 
respectively, the Fourier-transformed counterparts of CDM and baryon perturbations are given analytically as 

S l (k, t) = -3 [wxQeRi) - kR t cos(fci? 4 )] S z {t), (14) 

where i —c or b, and S c ^(t) are the mean overdensities of CDM and baryon within their respective top- hat regions. 
We can compute the neutrino perturbation in Fourier space from the zero-th moment of the perturbed distribution 
function: 

S v (k, r) = Ana^ 4 J q 2 dqef ^o{k, r). (15) 

The corresponding linear perturbation is given as S^(k, r) = Ana^ 4 J q 2 dqef^]^(k 1 r), a standard output of the CAMB 
code. 

We also need to compute the real-space overdensity of neutrino perturbations at each time step, which is needed 
for the spherical collapse model of CDM and baryon overdensities. The radial profile of neutrino density perturbation 
and its mass fluctuation within the CDM or baryon top-hat regions are 

<Ur;r) = f°° A7rk 2 dk^j^8„(k,T) , (16) 







Rc, b 

2 



6M V (< R c . b ;r) = / Anr A drp v 8 v {r), (17) 
Jo 

We use the publicly available code FFTLog [5l[ to compute the density profile at each time step, because the code 
allows for a fast computation of the integration involving the Bessel function kernel. Similarly, the density profile and 
the mass fluctuation for the linear neutrino perturbation can be obtained by using the linear density perturbation 
5^(fc,r) instead of 5^(k,T). 

We employ the decoupling epoch zjoc to set up the initial conditions of neutrino perturbations as in the baryon 
perturbations. Since the neutrino perturbations at Zdcc are well in the linear regime, we can determine the initial 
conditions by matching to the linear perturbation theory predictions. To be more precise, assuming the adiabatic 
initial conditions, we can determine the neutrino density perturbation at Zdcc around the top-hat CDM overdensity 
region [66j : 

Uk,T dec ) = Si(k,T dcc ) = l"l k ' TdCC U Uk,T dcc ), (18) 

where <5 c (fc, T dec ) is the Fourier transform of the CDM top-hat overdensity (Eq. P3]), and the functions T v (k, T dcc ) and 
T c {k-, idee) are the transfer functions of massive neutrinos and CDM at Zdcc, respectively. We use the CAMB outputs 
to obtain the transfer functions. The ratio T v /T c takes into account the relative amplitude difference of neutrino and 
CDM perturbations at Zdcc under the adiabatic initial conditions. We compute the zero-th moment of the perturbed 
distribution ^o{k, £dec) at the initial time by multiplying the CAMB output Tr[^Q (k, r dec )] with S c (k, Zdcc) /T c {k 1 T dcc ) 
so that Eq. (Til)]) gives the neutrino density perturbation around the CDM top-hat overdensity, where Tr[^>Q (k, T dcc )] 
is the linear transfer of the zero-th moment of the perturbed distribution function. Similarly, we compute the higher- 
order function ^i(k, Zdcc) (I > 1) from the CAMB outputs Tr[^^(A:, Zdec)] multiplied by S c (k 7 z dc c)/T c (k, z dcc ). We 
can obtain the subsequent evolution of neutrino perturbations by solving the Boltzmann equation hierarchies Eq. (|ll[) 
given these initial conditions at z dcc . 
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superhorizon 


subhorizon before the decoupling (zdec) 


subhorizon after Zdec 


CDM 
baryon 
massive- v 
radiation 


perturbed FRW (Eq. |SJ| <= 
perturbed FRW (Eq. [5J 
perturbed FRW (Eq[5j) 
perturbed FRW (Eq© 


=>■ sph. collapse model (Eq. [7|) 
5 b = 

5 U = 4 


sph. collapse model (Eq.O 
sph. collapse model (Eq. 1101) 
=> linearized Boltzmann eas. (Sec.lIIB3() 



TABLE I: A quick summary of our recipe used for solving the spherical top-hat collapse model in a multi-component system 
(CDM, baryon, massive neutrinos and radiation), where the initial top-hat CDM overdensity drives the collapse at late times. 
We set the initial time to be in the deeply radiation dominated regime, z% = 10 7 , in order to achieve a sufficient accuracy of 
setting up the initial conditions needed for the nonlinear dynamics. Also note that we assumed the adiabatic initial conditions 
which determine the relations between density perturbation amplitudes of different components at the initial time (see text 
for details). The density perturbations of interest are on superhorizon scales at the initial time, enter into the horizon and 
evolve on subhorizon scales. The notation denotes the matching of perturbations between different equations. Massive 

neutrinos with mass scales ( ^ a few 0.1 eV) become non-relativistic after the decoupling epoch (zd ec ), and we set up the 
neutrino perturbation amplitudes at z^ cc by matching with the CAMB outputs (see text for details). 

Our approach above is still an approximation; we used the linearized Boltzmann equations where each term in the 
hierarchies depends linearly on perturbation quantities. In other words, even though we include the effect of nonlinear 
gravitational potential, we ignore nonlinear terms in the Boltzmann equations, e.g., the term proportional to 0{^<j)), 
which can be important if the neutrino perturbation itself becomes nonlinear. We will come back to this issue later. 

C. Summary: our recipe for solving the spherical top-hat collapse model 

Here is a quick summary of the procedures we take for solving the spherical top-hat collapse model in a multi- 
component system of CDM, baryon and neutrinos. 

1. Choose a target mass scale of halo, M c , to determine the comoving scale of the spherical top-hat CDM overdensity 
region via the relation R C ^(M C ) = (3Af c /47rp C) o) 1 ^ 3 - 

2. Solve the linear CDM perturbation of the comoving scale i? Cj o from the initial time zi — 10 7 to the present 
time based on the linear perturbation theory, assuming the adiabatic initial conditions. In these calculations we 
properly take into account the fact that the density perturbations are on super-horizon scales in early epochs, 
enter into the horizon, and evolve on sub-horizon scales. 

3. Choose a target collapse redshift z co n that corresponds to the halo formation. Then, as a first guess, normalize 
the initial top-hat CDM overdensity, <5 c (z.;), for a given cosmological model in such a way that the linear-theory 
extrapolated overdensity satisfies the condition 8\(z co \\) = 1.686(1 + z C oii), a prediction for the collapse redshift 
for Einstein de-Sitter model, a CDM model without baryon and massive neutrino contributions. 

4. Solve the spherical top-hat collapse of CDM overdensity (Eq. [7]), coupled with the spherical top-hat collapse of 
baryon overdensity (Eq. |10| ) and the linearized Boltzmann equations of neutrino perturbations (see Sec. lIIB~3|) . 

5. Solve the nonlinear evolution of CDM overdensity until 8 C — > oo. Iteratively solve the spherical collapse model 
by changing the initial overdensity amplitude 8 c (zi) until the CDM overdensity becomes to collapse at the target 
collapse redshift, 8 c {z co \\) — > oo. Also obtain the linear-theory extrapolated overdensity for CDM or CDM plus 
baryon perturbations, a£ (z co n) or S^ b (z co \\). 

Table [fl also gives a quick summary of these treatments, clarifying which equations we use for solving the spherical 
top-hat collapse model. 

It would also be useful to explicitly list the assumptions we employ for the spherical collapse model: 

• We used the linearized Boltzmann hierarchy equations to solve the time evolution of linearized neutrino per- 
turbations. However, we include the effect of nonlinear gravitational potential due to the nonlinear CDM and 
baryon density perturbations. 

• We assumed the top-hat profiles of CDM and baryon perturbations. 

• We set the baryon and neutrino density perturbations to zero, i.e. 8^ = 8 U = 0, before the decoupling epoch, 
because we found it gives a good approximation compared to a more rigorous calculation under the adiabatic 
initial conditions. 
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For the second assumption above, rigorously speaking, even if we consider a spherically symmetric top-hat overdensity 
region at the initial time, the radial profile become changed by the presence of radiation, baryon and neutrino 
perturbations, which go out of the CDM overdensity region after the horizon crossing. As a result the overdensity 
region no longer obeys a top-hat profile. However, a spherical top-hat collapse model is anyway an approximated 
method for studying the nonlinear dynamics of the initial overdensity regions, preferentially representing density 
peaks in the primordial perturbations. Hence the top-hat overdensity region can be interpreted as the average density 
contrast around such density peaks after making a top-hat filtering with the smoothing scale of target halo mass 
scale [e.g.,l33l, for a similar discussion]. For these reasons, our treatment of assuming a spherically symmetric, top-hat 
overdensity for CDM perturbation is adequate enough for our purpose. Our main goal is to study the impact of 
massive neutrinos on the spherical collapse by comparing the results with and without massive neutrino contribution. 
Also note that our method includes the limit of spherical collapse for a pure CDM model without baryon and massive 
neutrino contributions, by imposing fib = Q v = 0. 



To compute the spherical collapse model, we assume a fl at-g eometry cold dark matter, A dominated cosmological 
model (ACDM) that is consistent with the WMAP results [43]. We further need to specify neutrino parameters. In 
this paper, we assume standard three flavors of neutrinos. Since structure formation is sensitive only to the sum 
of the three-species neutrino masses, we assume, for simplicity, that only one species of neutrinos are massive and 
other two species are massless. In this case, the neutrino free-streaming scale is shortest for a fixed total neutrino 
mass (or a fixed fi^o), and therefore the neutrino has the largest ability to cluster on small scales compared to a case 
that the total neutrino mass is split into different species. We then study how the spherical collapse is affected by 
massive neutrino assuming the mass scale ranging from 0.05 to a few 0.1 eV. This range of mass scales covers the 
lower and u ppe r bounds on the neutrino mass implied from the neutrino oscillation experiments and the cosmological 
constraints [l4j-[l6J. In most parts of this paper, when we add the massive neutrinos, we keep the energy density of 
"dark matter" (tt c o + f2„o) and other parameters fixed, and vary f2 c o. 

Fig. [T] shows how the CDM top-hat overdensity grows as a function of cosmic time. We assumed m„ = 0.05 eV for 
neutrino mass, a mass scale close to the lower bound implied from the normal mass hierarchy, and halo mass scale 
M = (47r/3)iZg ipej — 1O 14 /i _1 M0 corresponding to the comoving radius, i? Ci o = 6.89 /j _1 Mpc, for our fiducial ACDM 
model. By considering massive halos, we can estimate the largest effect of neutrinos on the spherical top-hat collapse, 
as such a massive halo has the ability to trap neutrinos around it due to the deepest gravitational potential well. We, 
as a working example, set the initial CDM overdensity so that the top-hat region collapses at redshift z co ii — 0.5 for 
a model without massive neutrino. 

The plot also shows the baryon density contrast within the CDM top-hat region as well as the radial profile of 
neutrino perturbations out to radii outside the top-hat region. Note that, for the baryon perturbation, we show only 
its density contrast within the CDM top-hat region (more exactly speaking, the baryon overdensity region is slightly 
more spatially-extended than the CDM as we described above). At sufficiently high redshifts such as z = 570, the 
baryon and neutrino perturbations are smaller in the amplitudes than the CDM density contrast. Then at lower 
redshifts, the baryon perturbation eventually catches up with the CDM perturbation. For this particular case, at 
lower redshifts z ?S 30, the comoving radius of CDM top-hat region starts to shrink, entering into the nonlinear regime 
(or equivalently deviating from the linear evolution). The figure explicitly demonstrates that the CDM and baryon 
perturbations, i.e. cold components, can collapse together having 5 C , St, — > oo at the collapse redshift. 

On the other hand, the neutrinos of this mass scale cannot catch up with the CDM perturbation due to the 
large free-streaming velocity. To be more precise, the present-day free-streaming scale in the comoving scale unit is 
Af s ~ (j v , v Hq 1 ~ 40 h~ 1 Mpc (see Appendix A in [HI), which is much larger than a few Mpc, a scale of the virial 
radius of massive halos. The plot shows that the neutrinos are indeed clustering around the CDM top-hat region, 
and become to have the radial profile. The neutrino perturbation peaks at the center of the CDM top-hat region, but 
the density contrast is still smaller than unity, so not yet in the highly nonlinear regime. More precisely, the neutrino 
perturbation averaged within the CDM top-hat region 8 U ~ 0.19 at the collapse redshift. We have also checked that, 
when neutrino mass is in the range smaller than a few 0.1 eV, the neutrino density contrast grows only up to the 
weakly nonlinear regime 8 V ~ a few even at the collapse time (for the case of m„,tot = 0.1 eV, 8 V ~ 0.51 at the collapse 
redshift). Note that, for the halo of 1O 15 M and the collapse redshift z co \\ ~ 0.5, 8 U ~ 0.54 and 2 for m^tot = 0.05 and 
0.1 eV, respectively). Therefore our approach using the linearized Boltzmann equations for neutrino perturbations is 
approximately validated. Once the halo is formed via virialization of the kinetic and gravitational bound energies, 
the neutrino would become stably clustered around the halo region as studied in [52J, [53j . Nevertheless the resulting 
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A. Spherical collapse in a mixed dark matter model 
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FIG. 1: Radial profile of density perturbation for each component (CDM, baryon and massive neutrino) in the CDM top-hat 
overdensity region assuming the adiabatic initial conditions to determine relative amplitudes of different components. We plot 
the profiles in units of comoving scale, Mpc/ft. We assumed m„ = 0.05 eV for neutrino mass scale, R = 6.89 /i~ 1 Mpc for radius of 
the top-hat region, which corresponds to halo mass scale M = 10 14 /i _1 M Q for our fiducial cosmological model, and determined 
the initial density amplitude so that the top-hat region collapses around 2 = 0.5. Note that, for the baryon perturbation, we 
show its density contrast within the CDM top-hat region for illustrative clarity (the baryon top-hat perturbation computed is 
spatially more extended, and holds the mass conservation within its own top-hat region). The different panels show the profiles 
at different redshifts as indicated. At sufficiently early redshift such as z = 570, after the decoupling epoch, <5b, 8 V <C <5 C , 
because baryon was coupled with radiation until the decoupling epoch (2 ~ 1100) and neutrino was free-streaming out of CDM 
potential well. As time goes by, the baryon perturbation eventually catches up with the CDM perturbation. Then at redshifts 
lower than z ~ 30 for this case, the top-hat radius starts to shrink and the top-hat dynamics deviates from the linear theory 
and enters into the nonlinear regime. As for neutrinos, the large velocity dispersion of neutrino particles prevents from catching 
up with the nonlinear collapse and then becomes to have a more spatially-extended profile than the CDM and baryon top-hat 
region. Even when CDM and baryon collapse (z ~ 0.5), the neutrino perturbation stays in the quasi nonlinear regime as <5„ < 1. 



neutrino overdensity is much smaller than the CDM and baryon perturbations in a halo region, therefore we ignore 
the neutrino mass contribution to the halo mass in the following analysis for simplicity. 

Fig- HI shows the time evolution of density contrasts within the top-hat region for each component. Note that we 
computed the neutrino density contrast by averaging the density profile within the top-hat region. The CDM and 
baryon collapse at z co n — 0.48 having S c , 5b — > 00. The neutrinos are affected by the nonlinear clustering of CDM 
perturbations, but do not enter into the highly nonlinear regime. 

In Fig. [31 we compare the spherical collapses of CDM perturbation for models with and without baryon perturbation 
and/or massive neutrino contributions. We used the same initial conditions of CDM perturbation except for the result 
without baryon perturbation (labelled as "w/o 5b and m v "). For the model without baryon perturbation (thin solid 
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FIG. 2: Time evolution of the density contrast of each component, averaged within the CDM top-hat region. We used the 
same initial conditions and model parameters as in Fig. [T] For this particular case, the neutrino perturbation averaged within 
the CDM top-hat region <5„ ~ 0.19 at the collapse redshift, and therefore is still in the quasi nonlinear regime. 

curve), we set the CDM perturbation amplitude to match the CDM amplitude at z& cc for our fiducial model (without 
massive neutrino), but set 51 m o = fi c o- In this case, if we set a sufficiently early collapse redshift, this model gives 
the collapse redshift given by # c (2 C oil) — 1-686(1 + z C oil)i the case for an Einstein de-Sitter model. At later collapse 
redshift, the cosmological constant becomes dominant in the cosmic expansion, and the collapse redshift differs from 
the Einstein de-Sitter prediction. This result is compared with other curves. First, the dashed curve shows the collapse 
of CDM perturbation i5 c when the baryon perturbation is added. The presence of baryon perturbation, which cannot 
grow during epochs from the horizon enter until the decoupling epoch, delays the spherical collapse. The delay is 
rather significant, from z co \\ ~ 0.9 to 0.5, and therefore the result suggests that we need to carefully take into account 
the effect of baryon perturbation on the nonlinear structure formation in N-body simulations, especially when we set 
up the initial conditions (we will discuss this issue later) . The bold solid and dotted curves show the spherical collapse 
when massive neutrino is further added for mass scales of 0.05 and 0.1 eV, respectively. The collapse redshifts are 
further delayed as z co u — 0.48 and 0.45, respectively. These mass scales correspond to the lower bounds on total 
neutrino mass for the normal and inverted mass hierarchies. Thus, adding the smoother, massive components into 
the CDM perturbation progressively delays the spherical collapse. 

We then study the linear-theory extrapolated overdensity at the collapse redshift (we will often call it the critical 
overdensity hereafter). In Fig. 2] we show the critical overdensity of CDM plus baryon perturbation, S^ b (z co \i; R), 
smoothed with length scales R — 6.89 and 14.8 ft, _1 Mpc corresponding to halo mass scales M = 10 14 and 10 15 h~ 1 M Ql 
respectively (the subscript "c6" stands for CDM plus baryon). The overdensity 5^ b (z; r) can serve as a clock to infer the 
collapse redshift, because it can be easily computed once the initial power spectra of CDM and baryon perturbations, 
halo mass scale and cosmological model are specified. In an Einstein de-Sitter universe (ft c o = 1), which includes CDM 
alone, the critical overdensity can be derived analytically [3(| Hl| and is found to be S]f(z c ) = 1.686, independently 
of halo mass and collapse redshift. Even for a model with curvature or dark energy contribution, the critical density 
Sc(z c ) ~ 1.686 to a good approximation, independently of halo mass (35rl37| . The top solid curve shows the critical 
overdensity when baryon is included, but massive neutrino is ignored. The critical overdensity differs from 1.686, 
and the change of ^(^colij R) is due to the presence of baryon perturbation, which has a smaller am plit ude than the 
CDM perturbation at higher redshifts (e.g., see Fig. [T]). Our result is consistent with the result in [34]], where they 
studied the effect of baryon perturbation on the spherical collapse at high redshift for much smaller halos that are 
relevant for first stars. Although the presence of smoother baryon perturbation delays the spherical collapse, it leads 
to the smaller S^ b (z co n; R) than 1.686. However, note that the linearly-extrapolated overdensity for CDM perturbation 
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FIG. 3: Time evolution of the CDM top-hat overdensity for different models. The dashed curve shows the result for our fiducial 
cosmological model without massive neutrino. We again assumed the halo mass scale M = 10 14 Mq and determined the 
initial CDM perturbation so that the top-hat region collapses at z ~ 0.5. The solid curve shows the result when ignoring the 
baryon perturbation, where the CDM perturbation amplitude is set so as to match that of the dashed curve at the decoupling 
epoch 2 — 1100. Note that the model of the solid curve leads the linear-theory extrapolated overdensity to be <5 L = 1.686 at 
redshifts before dark energy domination in the cosmic expansion (the result shown here is affected by dark energy domination). 
Comparing the solid and dashed curves manifests that the presence of baryon perturbation, which has a smaller amplitude at 
earlier redshifts as implied from Figs. [1] and [2l delays the spherical collapse. The solid and dotted curves show the results when 
further including massive neutrino for a fixed dark matter density f2 c o +f2„o- The neutrino perturbation is smoother than that 
of CDM perturbation, and delays the spherical collapse. 



alone, 8^ (z co i\; R), is indeed greater than 1.686. That is, the initial CDM top-hat overdensity with greater amplitude 
than expected from 5^(z co \i;R) = 1.686 is needed so that it collapses at a given collapse redshift z C oil- The curve 
peaks around z co \\ ~ 2 (a co ii — 0.33) having S^ b (z co \\; R) ~ 1.682. At the lower redshifts than z co \\ ~ 2, especially at 
z C o\\ ~ 1) the critical overdensity becomes smaller than the peak value due to the effect of the cosmological constant. 
When the cosmological constant or more generally dark energy becomes to dominate the cosmic energy density, the 
accelerating cosmic expansion slows down the growth of CDM plus baryon perturbation, and delays the spherical 
collapse. This yields the smaller critical overdensity. Both the linear and nonlinear growths of CDM plus baryon 
perturbation are delayed by the cosmic acceleration, and the linear growth is more suppressed than the nonlinear 
growth, because the spherical collapse eventually separates from the cosmic expansion in the nonlinear stage, and 
becomes more affected by the self-gravity of nonlinear CDM plus baryon overdensity. For these reasons, when the 
growth of density perturbations is suppressed by the faster cosmic expansion than the Einstein de-Sitter model, it 
generally leads to the smaller critical overdensity Sh(z co u; R) than 1.686. 

The other curves in Fig. @] show the results for 5^ b (z co \\; R) when including the massive neutrino for a fixed total 
matter density i7 m o- The presence of massive neutrino further delays the spherical collapse (see Fig. [3]), and in turn 
leads to the smaller critical overdensity S^ b (z co \\; R), the same trend for the effects of baryon perturbation and the 
cosmological constant. However, the massive neutrino only decreases S^ b (z co \\; R) by less than 0.1% compared to the 
solid curve, for these halo mass and neutrino mass scales and over a range of redshifts we have studied. This small 
change in S^ b (z co n; R) can be contrasted with the effect on the linear growth rate; the growth rate is suppressed by the 
amount of ~ 4/^ at relevant redshift compared to the growth rate without the massive neutrino (TTI . [l2[ , corresponding 
to 1.6 and 3.2% suppression for the neutrino mass scales of 0.05 and 0.1 eV, respectively. The results imply that the 
neutrino effect on the spherical collapse is well captured by the linear growth rate of CDM plus baryon perturbation. 

The different curves show that the change in 5^ b is not monotonic with changing neutrino masses, when keeping the 
present-day dark matter density (f2 c0 + fi„o) fixed. This non-trivial dependence can be understood as follows. The 
neutrino effect on the spherical collapse arises from its effect on the cosmic expansion history and the gravitational 
collapse of CDM perturbation. First, the presence of massive neutrino leads to a faster cosmic expansion during 
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FIG. 4: The linear-theory extrapolated density contrast of CDM plus baryon perturbation at the collapse redshift - the so-called 
critical density that can be used to infer the collapse redshift based on the linear theory. The left and right panels show the 
results for halo mass scales M = 10 14 and 10 15 /i _1 Mq, respectively. The different curves are the results without and with 
massive neutrino contribution assuming different neutrino mass scales. Note that the density contrast shown here is not for 
CDM perturbation alone, and the corresponding critical density of CDM perturbation is greater than shown in this plot. The 
overall change in the critical density from the Einstein de-Sitter result 8 L = 1.686 arises from the effect of baryon perturbation 
for higher redshifts, while the change at lower redshift z 1 is due to dark energy domination in the cosmic expansion. The 
effect of massive neutrino is in the range of the different curves. The curves show non-trivial dependence on neutrino mass 
scale (see the next figure). 



the neutrino was relativistic, which slows down the growth of CDM perturbation. Note that the neutrino becomes 
non-relativistic when Tomb — rn u . Secondly, the neutrino perturbation does contribute to the nonlinear gravitational 
collapse of CDM perturbation, and therefore accelerates the spherical collapse to some extent. The net effect arises 
from these competing effects. We can study which of these two effects is more important as follows. Fig. [5] shows 
the results when we ignore the neutrino perturbation (<$„ = 0) in both the spherical collapse calculation as well as 
the linear growth calculation for CDM plus baryon perturbation. The figure shows that, with increasing the neutrino 
mass scale, the spherical collapse more delays and the critical overdensity becomes monotonically smaller. Hence, 
comparing Figs. [4] and [5] manifests that the neutrino perturbation does contribute to the spherical collapse, which 
differs from the effect of smooth dark energy model [33, [54j • 

In Fig. O we summarize dependences of the critical overdensity on halo mass and neutrino mass, assuming the 
collapse redshift 2 co rj = 0. It can be found that, for a fixed halo mass scale, the critical overdensity first decreases with 
increasing the neutrino mass from 0.05 eV, but then starts to increase at greater neutrino mass scales from some mass 
scale. The turnover neutrino mass scale slightly changes with halo mass scale. This non-trivial dependence arises 
depending on which of the two competing effects discussed above dominates. If we ignore the neutrino perturbation, 
the critical overdensity decreases with increasing neutrino mass independently of halo mass. 

Summarizing the results in Figs. 21 [5] and [6j we can conclude that the effect of massive neutrino on the critical 
overdensity is very small, less than ~ 0.5%, compared to the critical overdensity without massive neutrino, for neutrino 
mass scales m u £ 0.5 eV and halo mass scales we are interested in. 



B. The impact of massive neutrinos on halo mass function 



In this subsection, we estimate the impact of massive neutrinos on the halo mass function that is one of the most 
important observables for cluster surveys. 

As we have shown, the effect of the massive neutrino on the nonlinear gravitational collapse of CDM plus baryon 
perturbation is well captured by the linear growth rate. In other words, the nonlinear neutrino clustering around the 
CDM overdensity does not largely change the nonlinear dynamics, and therefore is very unlikely to change structural 
properties of mass distribution within a halo. We here assume that the halo mass function for a MDM model can 
be obtained from a mapping of the mass function in CDM models without massive neutrino. That is, we assume 



13 




FIG. 5: Similarly to the previous figure, but we here ignored the effect of neutrino perturbation on the spherical collapse and 
on the linear density calculation; i.e. we set 8 V = and <5j; = 0. The critical overdensity becomes smaller at each redshift with 
increasing the neutrino mass scale. Therefore, comparing this figure with Fig. [4] clarifies that the non-trivial dependence of 5 oi , 
on neutrino mass scale is due to the effect of neutrino perturbation (see text for details). 



that the mapping of halo mass function can be obtained by assuming that (1) only the cold component (CDM plus 
baryon) can collapse to form halos, and (2) the halo mass function for a MDM model can be obtained just by replacing 
the linear-theory mass fluctuations appearing in the mass function for a CDM model with the corresponding mass 
fluctuation of CDM plus baryon perturbation for a MDM model: 
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(19) 



where the function f(v) is the fitting formula that is obtained based on a suit of N-body simulations for CDM 
models. The previous works have shown that the fitting formula is well characterized in terms of the peak height, 
v = 8^ 1 . it /cr(M; z) (UUll, where 6^ lit is the linear-theory extrapolated critical overdensity for halo formation at a 
given redshift and cr(M; z) is the linear rms mass fluctuation smoothed with the halo mass scale M and at redshift 
z. In Eq. (jT9"]) . we assumed that we can obtain the halo mass function for a MDM model simply by using the peak 
height for CDM plus baryon perturbation as well as by using the prefactor p c i,/M, the mean mass density of CDM 
plus baryon, because the cold component is the collapsing component to form halos. To be more precise, the rms 
mass fluctuation of halo mass M for CDM plus baryon perturbation is defined as 
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where P^ b (k; z) is the linear power spectrum of CDM plus baryon perturbation at target redshift z, and W(kR,M) the 
Fourier-transformed top-hat filter: W(x) = 3(sinir — a; cos a;) /a; 3 . The filtering scale and halo mass are related via 
M = {Air /3)p c bfiR\j (p c b,o is the present-day mean mass density of CDM plus baryon). 

As for the fitting formula, we use the formula in [lj| that is obtained from N-body simulations for a range of CDM 
models varying around the fiducial cosmological model consistent with the WMAP data: 



f(y,z) = A\j -exp 
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(21) 



where A = 0.333(1 + z)" ' 11 , a = 0.788(1 + z)~ om , p = 0.807 and q = 1.795. For the peak height v = 6 crit /a(M; 2), 
[Tsj simply used the fixed critical overdensity <5 cr it = 1.686, the value of Einstein de-Sitter model, and then found 
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FIG. 6: The plot shows how the critical overdensity 8^ b changes with different halo mass scales and neutrino mass scales. Here 
we consider z co ii = for the collapse redshift. For comparison, the dashed curve shows the result without massive neutrino. 
More massive halos have a deeper gravitational potential well and therefore are more capable of capturing neutrinos around it. 
In addition, neutrinos of greater mass scales have a smaller free-streaming scale, and are more captured by the CDM top-hat 
region. These facts together with Figs. 2] and \5\ explain the non-trivial dependence of 8^ b on halo mass scales as well as neutrino 
mass scales. 

the best-fit parameters A, a and so on by fitting the functional form above with the mass function measured from 
simulations for variant CDM models. More exactly speaking, when we have the effects of baryon perturbation and 
cosmic acceleration, the critical overdensity 5 h is changed from the Einstein de-Sitter value 5 L = 1.686. However, 
the change is very small, less than a percent level (see Fig. and therefore it was assumed that the change of S L is 
absorbed by tuning the fitting model parameters. If the change of <5 L is properly taken into account, the fitting will 
yield slightly different best-fit model parameters of A, a and so on. Furthermore, although the presence of bar yon 
perturbation changes the collapse of CDM perturbation (see Fig. [3]), we here assume that the simulations in [18[ 
properly take into account the effect of baryon when setting up the initial conditions of N-body simulations (see below 
for a further discussion). To compute a c b(M;z) in Eq. (flU|) , we use the CAMB code [5(| to compute the transfer 
functions. The CAMB outputs include the effect of massive neutrinos or baryon perturbations on the growth of CDM 
perturbation. 

The upper panel of Fig. [7] compares the halo mass functions at z = for different models with and without massive 
neutrino contribution. The lower panel explicitly shows the ratio of the mass functions with and without massive 
neutrino. Here we assumed m u — 0.05 and 0.1 eV for the neutrino mass scale, which are close to the lower bounds 
of the normal mass hierarchy (NH) and the inverted hierarchy (IH) that are implied from the terrestrial experiments. 
Hence either of these results would inevitably exist in our universe. The presence of massive neutrino decreases the 
abundance of halos, more significantly for more massive halos that reside in the exponential tail of mass function. 
The decrease in the halo abundance is up to a factor 2 around ~ 5 x 10 15 }i~ 1 Mq, This change can be compared 
to the effect on the linear mass fluctuation such as as] the neutrino of these mass scales decreases as only by a few 
percent for neutrinos of these mass scales. Again the higher sensitivity of halo mass function to neutrino mass is 
through the exponential tail of mass function at massive halo ends. The thin solid curve in the lower panel (although 
almost overlapped with the dotted curve) shows the ratio when further taking into account the change in the critical 
density 8^ h in the mass function (Eq. [H?)) ; more explicitly we decreases the critical density by 0.03%, a maximum 
change implied from Fig. [5] for the case of m u xot — 0.1 eV. It is clear that the change in the critical density due to 
the massive neutrinos causes a negligible effect on the halo mass function. 

Fig.[5]shows the similar results, but for higher redshifts z = 1 and 1.5, respectively. The decrease in the abundance 
of cluster scale-halos is more significant at higher redshifts. 

One may think whether or not the effect of massive neutrino on the halo mass function is mostly described by the 
change of as, the normalization parameter of power spectrum amplitudes often used in the literature. Fig. [9] compares 
the halo mass functions for a MDM model with m v = 0.1 eV and for a ACDM model where as at z — 1 is lowered so 
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FIG. 7: The upper panel shows the halo mass function at z = for CDM-dominated models with and without massive 
neutrinos. The halo mass function for a mixed dark matter model (CDM plus massive neutrino) is computed by mapping the 
fitting formula for CDM model based simulations using Eq. (|19p . For the solid and dotted curves, we assume the neutrino 
mass scales m„ = 0.05 and 0.1 eV, which are close to the lower mass bounds for the normal and inverted mass hierarchies (NH 
and IH), respectively. The presence of massive neutrino, for a fixed £7 c q + n„o, decreases the abundance of massive halos. The 
lower panel explicitly shows the ratio of the mass functions for models with and without massive neutrino contribution. The 
linear mass fluctuation such as erg changes only by a few percent at most for these neutrino masses, however, the abundance of 
massive halos may decrease by up to a factor 2 at a few 10 15 /i -1 Mq. 




as to match the as value for the MDM model; more precisely, &s(z — 1) is changed to 0.486 from 0.500. Note that 
both the models have the same fi c o + 0„q. The ACDM model with the lowered as roughly reproduces the decrease 
in the halo abundance. However, the two curves do not exactly agree because of the difference in the linear power 
spectra of CDM and baryon perturbations. Nevertheless, this gives a justification of the neutrino mass constraint 
derived in [lj, where the neutrino mass constraint is obtained from the allowed range of as values that are derived 
by comparing the observed abundance of X-ray luminous clusters with the model halo mass functions varying within 
CDM models without massive neutrino contribution. 
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FIG. 9: Similar plot to the lower panel in the previous figure, but we here compare the halo mass function for a MDM model 
with m„ = 0.1 eV to the mass function for a CDM model (without massive neutrino), where as value is lowered by the same 
amount as in the neutrino suppression effect on the linear mass fluctuation at 8 /i _1 Mpc for the MDM model. More exactly, 
the a 8 value at z = 1 is changed to 0.486 from 0.500 for the dashed curve. The CDM model with normalization of the lowered 
ag value reproduces the mass function for the MDM model, for the same S7 c o + fi„o, within 30% level accuracy over the range 
of halo masses we consider. 
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FIG. 10: As in Fig. [4] the linear-theory extrapolated critical density for halos of M = 10 15 /i _1 A/q as a function of collapse 
redshift for MDM models of different neutrino mass scales, where we added massive neutrinos around the fiducial ACDM 
model by changing either h or Qa parameter with fixing the CDM density parameter fl c oh 2 and keeping the flat geometry 
fico + fibo + £2i/0 + = 1 (see text for details). Note that, for the previous plots, we varied the CDM dark matter density fi c o 
with fixing the total dark matter density Q c q + fi„o to the fiducial value when adding massive neutrinos. 



C. Discussion: Cosmological parameter degeneracies 

When adding the massive neutrinos for different masses, we have so far kept the total dark matter density, r2 c o + f2yOi 
fixed. For a more practical perspective, the CMB information give precise constraints on the CDM and baryon 
densities, ri c o^ 2 and flhoh 2 , as well as the curvature parameter or equivalently the total energy density, Sl c o + f2bo + 
f2„o + Qa — 1. Massive neutrinos with small mass scales of a few 0.1 eV were relativistic before the decoupling 
epoch, and do not affect the CMB observables. Therefore the CMB observables cannot well constrain the neutrino 
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FIG. 11: Shown is the ratio of the halo mass functions with and without massive neutrinos of m„ jt ot = 0.1 eV. The different 
curves show the results when changing either h, J1a or Q c o alone around our fiducial cosmological model (see the end of Sec. HI 
with fixing the parameters flboh 2 and assuming a flat geometry of f2 c o + fibo + QvO + = 1. We have so far considered 
the case varying Q c o, and the cases varying h or Qa are motivated by the fact that the CMB observables well constrain the 
curvature parameter, il c h 2 and flboh 2 . Around our fiducial cosmological model, these are equivalent to the parameter changes, 
h = 0.71 ->• h = 0.7128, tt A = 0.7354 ->■ 0.7333 or Q c0 = 0.2198 ™s> 0.2177, respectively. 

mass of the small mass scales, leaving degeneracies in cosmological parameters. Given these facts one might think 
that, when adding the massive neutrinos, we should keep these CMB-constrained parameters fixed. If we assume a 
fiat geometry, this is equivalent to varying either the Hubble parameter h or the energy density of the cosmological 
constant 51a with fixing the CMB parameters above. For example, when the neutrino mass m^tot =0.1 eV is added, 
this leads to h = 0.7128 or = 0.7333 from the fiducial values h — 0.71 or ^a = 0.7354, respectively. 

Fig. HOl shows the critical density for halos of M — 15 15 ft.~ 1 M Q for a MDM model with various neutrino mass scales, 
where we varied either Q,^ or h by the amount determined by the neutrino mass scale, but fixing Q. c0 h 2 . The results 
are similar to Fig. @] the effect of massive neutrinos on the critical density is very small for the range of cosmological 
models. Fig. [TT] shows how the MDM models alter the halo mass function compared to the case without massive 
neutrinos, which can be compared with our fiducial case where the massive neutrinos of m v tot =0.1 eV are added 
by varying the CDM density parameter S7 c o - The parameter change of h or ft a also leads to the smaller abundance 
of massive halos as in the case changing f2 c o, but the decrease is slightly smaller than the case when changing O c o- 

IV. SUMMARY 

In this paper, we have developed a method to solve the nonlinear dynamics of top-hat CDM overdensity region 
including the effects of baryon perturbation and massive neutrinos. In developing the spherical collapse model, 
we properly set up the initial conditions of each components (baryon, CDM and neutrinos), which have different 
amplitudes and profiles, assuming the adiabatic initial conditions (see Fig. [1}. In fact we found that the nonlinear 
dynamics is very sensitive to detailed setup of the initial conditions of top-hat CDM perturbation, more precisely 
5 c (zi) and the velocity of top-hat radius R{z{). For example, we cannot employ the linear-theory prediction for an 
Einstein de-Sitter model, 8 oc a, to set up the initial conditions, e.g., even at an epoch in the sufficiently linear regime 
such as the decoupling epoch z- m { ~ 1100, because this solution ignores that the CDM perturbation is affected by the 
presence of baryon and massive neutrino. 

Since we cannot treat the neutrinos as a perfect fluid, we properly solved the linearized Boltzmann hierarchy 
equations to compute time evolution of linearized neutrino perturbations, where we include the effect of nonlinear 
gravitational potential due to the nonlinear CDM and baryon perturbation in the late stage. For neutrino mass 
scales lighter than a few 0.1 eV, the range inferred from the neutrino oscillation experiments and the cosmological 
constraints, the neutrino perturbation stays in the quasi nonlinear regime, 8 y ^ 1 (see Fig.[^J). This gives & justification 
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of our treatment where we used the linearized Boltzmann equations. As for an improved modeling, one can further 
include the nonlinear terms such as the coupling term between the nonlinear gravitational potential and the perturbed 
phase-space density of neutrinos in order to solve the time evolution of neutrino perturbations in a perturbation theory 
manner. 

By solving the spherical collapse model for cosmological models around a ACDM model that is consistent with 
the WMAP data, we found that both the neutrino and baryon perturbations delay the collapse of CDM overdensity 
compared to a model with CDM alone (Fig. However, interestingly we found that the collapse redshift can be 
well monitored by the linear-theory extrapolated overdensity of CDM (plus baryon) perturbation(s) for the ranges of 
neutrino masses ( a few 0.1 eV) and halo mass scales we have considered. This result is promising because the 
linear-theory extrapolated overdensity (the critical density) can be accurately computed using the linear perturbation 
theory, once cosmological model and neutrino mass are specified. In other words, we found that the massive neutrinos 
with the rang of mass scales lead to only a small change in the critical density by < 0.1% compared to the model 
without massive neutrino, but with the same f2 m o (Figs. [4l [5j and|6|). 

Given the results of the spherical collapse model, we gave Eq. (fl~9l) to estimate the halo mass function including the 
effect of massive neutrinos, where the effect of massive neutrinos are properly taken into account in the linear mass 
fluctuations of CDM and baryon perturbations at a given redshift, smoothed with a given halo mass scale; cr c f,(M, z) 
[also seel23l for the similar discussion]. Using the equation, we found that the presence of massive neutrinos with 0.05 
and 0.1 eV, the lower-bound mass scales of normal and inverted mass hierarchies, respectively, may cause a significant 
decrease in the abundance of massive halos; more specifically, up to a factor of 2 for halos with 10 15 M Q and at z ~ 1 
(see Figs. [7] and IS]). Thus our results imply that massive neutrinos, which should exist in our universe, relax to some 
extent a possible tension that the cutting-edge SZ experiments could not find as many massive clusters as what was 
originally expected @, H| ■ This needs to be further studied more carefully. Since it is still challenging to accurately 
simulate nonlinear structure formation in a MDM model, especially for such light neutrino mass scales of < a few 
0.1 eV [see [221423, for the attempts], the analytical model developed in this paper will give a useful tool or at least 
useful guidance for interpreting ongoing and upcoming wide-area surveys of massive clusters. 

Our findings also propose several applications. First, as we stressed above, a careful setup of the initial conditions is 
very important in order to have an accurate nonlinear dynamics, for a multi-component system with CDM, baryon and 
neutrinos. This implies that it is very important to set up the accurate initial conditions for cosmological simulations 
including the effect of baryon such as smoothed particle hydrodynamical (SPH) simulations [see l56l458l . for the similar 
discussion]. Since the spherical collapse model gives an exact solution of the nonlinear dynamics, albeit an unrealistic 
symmetry assumed, we can explore how to set up the initial conditions by combining the spherical collapse model with 
the linear and/or perturbation theory predictions. For example, it was shown that using the second-order Lagrangian 
perturbation theory allows one to set up more accurate initial conditions of N-body simulations that are simulations 
for a model with CDM-alone or single cold component [59|, [6(| • We can extend this analysis to a multi-component 
system; we can apply the second-order Lagrangian perturbation theory to CDM and baryon perturbations separately 
by taking into account the different growth rates, and then can study how the improved initial conditions can reproduce 
the exact solution of spherical collapse for CDM and baryon perturbations starting from a given initial redshift. Such 
a study will give a useful guidance for exploring how to set up the initial conditions for CDM and baryon particles 
in a SPH-type simulation. This can be further extended to a case further including the neutrino particles. These are 
our future study, and will be presented elsewhere. 

Secondly, several studies recently claimed that detected massive clusters at high redshifts beyond z ~ 1 may give 
a tension of ACDM structure formation model [HJ, H2, also see references therein] . It is indeed interesting to explore 
whether or not these particular catalogs of clusters, which are found by different observations/surveys under different 
selection functions, can falsify the ACDM predictions as explored in [63|. However, the effect of massive neutrinos 
has been ignored in the previous studies. Again, the method developed in this paper can be used to address how the 
presence of the high-z massive clusters may falsify a more realistic cosmological model that includes massive neutrino 
contribution. This study will be presented elsewhere. 
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